# Importing the required libraries
import numpy as np
import scipy.stats as st
import scipy.signal as sg
import matplotlib.pyplot as plt
import scipy.ndimage as nd
import cv2
# To remove unnecessary warnings
import warnings
warnings.filterwarnings("ignore")
# Note : Used for autoformatting the notebook, remove if causing any error
%load_ext nb_black
# Parameters for good matplotlib output
plt.rcParams["figure.autolayout"] = False
plt.rcParams["figure.figsize"] = 9, 6
plt.rcParams["axes.labelsize"] = 18
plt.rcParams["axes.titlesize"] = 20
plt.rcParams["font.size"] = 16
plt.rcParams["lines.linewidth"] = 2.0
plt.rcParams["lines.markersize"] = 6
plt.rcParams["legend.fontsize"] = 18
plt.rcParams["legend.numpoints"] = 2
plt.rcParams["legend.loc"] = "best"
plt.rcParams["legend.fancybox"] = True
plt.rcParams["legend.shadow"] = True
plt.rcParams["text.usetex"] = True
plt.rcParams["font.family"] = "serif"
plt.rcParams["font.serif"] = "cm"
plt.rcParams["text.latex.preamble"] = [r"\usepackage{amsmath}", r"\usepackage{amssymb}"]
Creating the signal generator using numpy
# Function to generate signal X given range k
def x(k):
x = np.zeros(k.shape)
x[(k > -1) & (k < 16)] = 3 + np.sin(2 * k[(k > -1) & (k < 16)] * np.pi / 15)
return x
Plotting the signal
# Creating the signal range
k = np.arange(0, 16)
# Generating the signal and plotting it
plt.title(r"Input Signal $x_k$")
plt.plot(k, x(k), "bo")
plt.xlabel("k")
plt.ylabel("x")
plt.show()
def y_1(k):
return x(k + 1) - x(k)
plt.title(r"Filtered Signal $y_k = x_{k+1}-x_k$")
plt.plot(k, y_1(k), "bo")
plt.xlabel("k")
plt.ylabel("y")
plt.show()
def y_2(k):
return x(k) - np.average(x(k))
plt.title(r"Filtered Signal $y_k = x_{k}-\bar{X}$")
plt.plot(k, y_2(k), "bo")
plt.xlabel("k")
plt.ylabel("y")
plt.show()
def y_3(k):
return np.array([np.median(x(np.arange(l - 2, l + 3))) for l in k])
plt.title(r"Filtered Signal $y_k = median(\{x_l : l \in [k-2,k+2]\})$")
plt.plot(k, y_3(k), "bo")
plt.xlabel("k")
plt.ylabel("y")
plt.show()
Using linear interpolation :
$x_{k+0.5} = x_{k}+\dfrac{(x_{k+1} - x_{k})}{2} = \dfrac{x_{k+1}+x_{k}}{2} \implies$ The values at the middle (0.5) are just the averages of the nearby points.
$\implies y_k = \dfrac{x_{k}+x_{k+1} - x_{k}-x_{k-1}}{2} = \dfrac{x_{k+1}-x_{k-1}}{2}$
def y_4(k):
return (x(k + 1) - x(k - 1)) / 2
plt.title(r"Filtered Signal $y_k = x_{k+0.5}-x_{k-0.5}$")
plt.plot(k, y_4(k), "bo")
plt.xlabel("k")
plt.ylabel("y")
plt.show()
(We use the same transform as above)
def y_5(k):
return np.abs((x(k + 1) - x(k - 1)) / 2)
plt.title(r"Filtered Signal $y_k = |x_{k+0.5}-x_{k-0.5}|$")
plt.plot(k, y_5(k), "bo")
plt.xlabel("k")
plt.ylabel("y")
plt.show()
def y_6(k):
return np.array([np.sum(x(np.arange(l - 2, l + 3))) / 5 for l in k])
plt.title(r"Filtered Signal $y_k = \frac{1}{5}\sum_{i=k-2}^{k+2} x_i$")
plt.plot(k, y_6(k), "bo")
plt.xlabel("k")
plt.ylabel("y")
plt.show()
Defining function to do 1D convolution of two signals.
def convolve_1d(x, h):
# Switching if length of x is lesser than h
if x.shape[0] < h.shape[0]:
temp = x
x = h
h = temp
# Padding x for handling edge cases
x_padded = np.hstack([np.zeros((len(h) // 2,)), x, np.zeros((len(h) // 2,))])
# Reversing the filter for convolution operation
h = h[::-1]
y = []
# Convolution operation
for i in range(len(x)):
f = 0
for k in range(len(h)):
f += x_padded[i + k] * h[k]
y.append(f)
return np.array(y)
The convolutional filter for $y_1 \to h_1 = \begin{bmatrix}1 & -1 & 0\end{bmatrix}$
h_1 = np.array([1, -1, 0])
The convolutional filter for $y_4 \to h_4 = \begin{bmatrix}0.5 & 0 & -0.5\end{bmatrix}$
h_4 = np.array([0.5, 0, -0.5])
The convolutional filter for $y_6 \to h_6 = \begin{bmatrix}0.2 & 0.2 & 0.2 & 0.2 & 0.2\end{bmatrix}$
h_6 = np.array([0.2, 0.2, 0.2, 0.2, 0.2])
Applying the filter $h_1$ and visualizing the output and comparing with direct computation
y_1_c = convolve_1d(x(k), h_1)
fig, axs = plt.subplots(1, 2, figsize=(14, 6))
plt.suptitle(r"$y_k = x_{k+1}-x_k$")
axs[0].plot(k, y_1(k), "bo")
axs[0].set_title("Original")
axs[0].set_xlabel("k")
axs[0].set_ylabel("y")
axs[1].plot(k, y_1_c, "bo")
axs[1].set_title("Using Convolution")
axs[1].set_xlabel("k")
axs[1].set_ylabel("y")
plt.show()
Applying the filter $h_4$ and visualizing the output and comparing with direct computation
y_4_c = convolve_1d(x(k), h_4)
fig, axs = plt.subplots(1, 2, figsize=(14, 6))
plt.suptitle(r"$y_k = x_{k+0.5}-x_{k-0.5}$")
axs[0].plot(k, y_4(k), "bo")
axs[0].set_title("Original")
axs[0].set_xlabel("k")
axs[0].set_ylabel("y")
axs[1].plot(k, y_4_c, "bo")
axs[1].set_title("Using Convolution")
axs[1].set_xlabel("k")
axs[1].set_ylabel("y")
plt.show()
Applying the filter $h_6$ and visualizing the output and comparing with direct computation
y_6_c = convolve_1d(x(k), h_6)
fig, axs = plt.subplots(1, 2, figsize=(14, 6))
plt.suptitle(r"$y_k = \frac{1}{5}\sum\limits_{i=k-2}^{k+2} x_i$")
axs[0].plot(k, y_6(k), "bo")
axs[0].set_title("Original")
axs[0].set_xlabel("k")
axs[0].set_ylabel("y")
axs[1].plot(k, y_6_c, "bo")
axs[1].set_title("Using Convolution")
axs[1].set_xlabel("k")
axs[1].set_ylabel("y")
plt.show()
We will try to transform the convolution in space domain into multiplication in frequency domain and transform back to time domain to check the output.
Computing the convolution using multiplication in frequency domain for $h_1$ and visualizing the output
y_1_ft = np.fft.ifft(
np.fft.fft(np.concatenate([h_1, np.zeros(len(k) - len(h_1))])) * np.fft.fft(x(k))
)
fig, axs = plt.subplots(1, 2, figsize=(14, 6))
plt.suptitle(r"$y_k = x_{k+1}-x_k$")
axs[0].plot(k, y_1(k), "bo")
axs[0].set_title("Original")
axs[0].set_xlabel("k")
axs[0].set_ylabel("y")
axs[1].plot(k, y_1_ft, "bo")
axs[1].set_title("Using Fourier Transoform")
axs[1].set_xlabel("k")
axs[1].set_ylabel("y")
plt.show()
Computing the convolution using multiplication in frequency domain for $h_4$ and visualizing the output
y_4_ft = np.fft.ifft(
np.fft.fft(np.concatenate([h_4, np.zeros(len(k) - len(h_4))])) * np.fft.fft(x(k))
)
fig, axs = plt.subplots(1, 2, figsize=(14, 6))
plt.suptitle(r"$y_k = x_{k+0.5}-x_{k-0.5}$")
axs[0].plot(k, y_4(k), "bo")
axs[0].set_title("Original")
axs[0].set_xlabel("k")
axs[0].set_ylabel("y")
axs[1].plot(k, y_4_ft, "bo")
axs[1].set_title("Using Fourier Transoform")
axs[1].set_xlabel("k")
axs[1].set_ylabel("y")
plt.show()
Computing the convolution using multiplication in frequency domain for $h_6$ and visualizing the output
y_6_ft = np.fft.ifft(
np.fft.fft(np.concatenate([h_6, np.zeros(len(k) - len(h_6))])) * np.fft.fft(x(k))
)
fig, axs = plt.subplots(1, 2, figsize=(14, 6))
plt.suptitle(r"$y_k = \frac{1}{5}\sum\limits_{i=k-2}^{k+2} x_i$")
axs[0].plot(k, y_6(k), "bo")
axs[0].set_title("Original")
axs[0].set_xlabel("k")
axs[0].set_ylabel("y")
axs[1].plot(k, y_6_ft, "bo")
axs[1].set_title("Using Fourier Transoform")
axs[1].set_xlabel("k")
axs[1].set_ylabel("y")
plt.show()
We see that none of the above outputs from Fourier Transform match the output of applying the filter. This is because the multiplication in frequency domain does not match directly the linear convolution in space domain, rather it matches circular convolution in space domain. To get normal linear convolution from circular convolution, we need to first pad the input signal by half the length of the filter signal on both sides, then multiply the Fourier Transforms of the signals and then clip out the starting part of the output which is not required to get the required output.
# Function to do linear convolution using fourier transform outputs
def fft_convolve_1d(x, h):
x_padded = np.hstack([np.zeros(len(h) // 2,), x, np.zeros(len(h) // 2,)])
y = np.fft.ifft(
np.fft.fft(np.concatenate([h, np.zeros(len(x_padded) - len(h))]))
* np.fft.fft(x_padded)
)[(len(h) // 2) * 2 :]
return y
Applying the modified convolution using Fourier Transform and visualizing the output for $h_1$
y_1_fft_con = fft_convolve_1d(x(k), h_1)
fig, axs = plt.subplots(1, 2, figsize=(14, 6))
plt.suptitle(r"$y_k = x_{k+1}-x_k$")
axs[0].plot(k, y_1(k), "bo")
axs[0].set_title("Original")
axs[0].set_xlabel("k")
axs[0].set_ylabel("y")
axs[1].plot(k, y_1_fft_con, "bo")
axs[1].set_title("Using Fourier Transoform(Changed)")
axs[1].set_xlabel("k")
axs[1].set_ylabel("y")
plt.show()
Applying the modified convolution using Fourier Transform and visualizing the output for $h_4$
y_4_fft_con = fft_convolve_1d(x(k), h_4)
fig, axs = plt.subplots(1, 2, figsize=(14, 6))
plt.suptitle(r"$y_k = x_{k+0.5}-x_{k-0.5}$")
axs[0].plot(k, y_4(k), "bo")
axs[0].set_title("Original")
axs[0].set_xlabel("k")
axs[0].set_ylabel("y")
axs[1].plot(k, y_4_fft_con, "bo")
axs[1].set_title("Using Fourier Transoform(Changed)")
axs[1].set_xlabel("k")
axs[1].set_ylabel("y")
plt.show()
Applying the modified convolution using Fourier Transform and visualizing the output for $h_6$
y_6_fft_con = fft_convolve_1d(x(k), h_6)
fig, axs = plt.subplots(1, 2, figsize=(14, 6))
plt.suptitle(r"$y_k = \frac{1}{5}\sum\limits_{i=k-2}^{k+2} x_i$")
axs[0].plot(k, y_6(k), "bo")
axs[0].set_title("Original")
axs[0].set_xlabel("k")
axs[0].set_ylabel("y")
axs[1].plot(k, y_6_fft_con, "bo")
axs[1].set_title("Using Fourier Transoform(Changed)")
axs[1].set_xlabel("k")
axs[1].set_ylabel("y")
plt.show()
Importing the function from helpers.py
from helpers import load_image, save_image, vis_hybrid_image
Finding the best method to do 2D convolution for given question
img1 = load_image("data/ex01/bicycle.bmp")
sg.choose_conv_method(img1, np.ones((3, 3)))
Since the output above is given as 'fft', we will use fft for convolution in 2D case
Function to do 2D convolution using FFT
def convolve_2d(img, fil):
# Finding the padding value for the filter
## Horizontal Padding
hpad1 = (img.shape[0] - fil.shape[0]) // 2
## Handing even/odd cases
if (img.shape[0] - fil.shape[0]) % 2:
hpad2 = hpad1 + 1
else:
hpad2 = hpad1
## Vertical Padding
vpad1 = (img.shape[1] - fil.shape[1]) // 2
## Handling even/odd cases
if (img.shape[1] - fil.shape[1]) % 2:
vpad2 = vpad1 + 1
else:
vpad2 = vpad1
# Padding the filter
fil_new = np.pad(fil, [(hpad1, hpad2), (vpad1, vpad2)],)
# Transforming the image and filter to Fourier Domain
fr1 = np.fft.fft2(img)
fr2 = np.fft.fft2(
np.flipud(np.fliplr(fil_new))
) # Flipping the filter to conform to convolution similar to CNN
m, n = fr1.shape
# Finding the convolution output in Fourier Domain
cc = np.real(np.fft.ifft2(fr1 * fr2))
# Transforming the output to match the actual convolution output
cc = np.roll(cc, -m // 2 + 1, axis=0)
cc = np.roll(cc, -n // 2 + 1, axis=1)
# Clipping the outputs
cc[cc > 1] = 1
cc[cc < 0] = 0
return cc
Function to convolve an image with a given filter
def filter_image(img, fil):
# Checking for correctness of filter
if fil.shape[0] % 2 == 0 or fil.shape[1] % 2 == 0:
raise Exception("Please enter correct filter")
# Applying filter on color image using the convolution function defined above
if len(img.shape) != 2 and img.shape[2] == 3:
out = []
for i in range(3):
out.append(convolve_2d(img[:, :, i], fil))
return np.dstack(out)
# Appying filter on grayscale image using the convolution function defined above
else:
return convolve_2d(img, fil)
Function to create a gaussian kernel
def gaussian_kernel(sig):
# Creating kernel size for given sigma value
if (int(6 * sig) + 1) % 2 == 0:
kern = int(6 * sig)
else:
kern = int(6 * sig) + 1
# Creating an 1d Space
x = np.linspace(-2 / sig, 2 / sig, kern)
# Generating 1d Gaussian Kernel by taking the difference of CDF of the above space
gauss_1d = st.norm.pdf(x)
# Creating the 2d Gaussian Kernel by multiplying two 1d Gaussian Kernels
gauss_2d = np.outer(gauss_1d, gauss_1d)
# Normalizing the Gaussian Kernel
gauss_kernel = gauss_2d / gauss_2d.sum()
return gauss_kernel
We will define functions to create the hybrid filters :
low_pass : Gives the low-pass filtered output of image.high_pass : Gives the high-pass filtered output of image.hybrid_image : Combines the above two functions to create the hybrid image.def low_pass(img, freq):
# Defining the Gaussian kernel
kernel = gaussian_kernel(freq)
# Filtering the image with above kernel
out = filter_image(img, kernel)
# Clipping and normalizing the output
out = out / out.max()
out[out < 0] = 0
out[out > 1] = 1
return out
def high_pass(img, freq):
# Defining the Gaussian kernel
kernel = gaussian_kernel(freq)
# Filtering the image with above kernel
out = img - filter_image(img, kernel)
# Clipping and normalizing the output
out = out / out.max() + 0.5
out[out < 0] = 0
out[out > 1] = 1
return out
def hybrid_image(img1, img2, freq1, freq2):
# Finding the low_pass and high_pass images
img_low = low_pass(img1, freq1)
img_high = high_pass(img2, freq2)
# Combining the high_pass and low_pass
return (img_low + img_high) / 2
Bicycle and Motorcycle
img1 = load_image("data/ex01/bicycle.bmp")
img2 = load_image("data/ex01/motorcycle.bmp")
plt.figure(figsize=(15, 10))
plt.title("Bicycle and Motorcycle")
plt.imshow(np.hstack([img1, img2]))
plt.axis("off")
plt.show()
img_hyb = hybrid_image(img1, img2, 5, 3)
plt.figure(figsize=(15, 10))
plt.title("Highpass Motorcycle and Lowpass Bicycle")
plt.imshow(vis_hybrid_image(img_hyb))
plt.axis("off")
plt.show()
img_hyb = hybrid_image(img2, img1, 3, 1.5)
plt.figure(figsize=(15, 10))
plt.title("Highpass Bicycle and Lowpass Motorcycle")
plt.imshow(vis_hybrid_image(img_hyb))
plt.axis("off")
plt.show()
Bird and Plane
img1 = load_image("data/ex02/bird.bmp")
img2 = load_image("data/ex02/plane.bmp")
plt.figure(figsize=(15, 10))
plt.title("Bird and Plane")
plt.imshow(np.hstack([img1, img2]))
plt.axis("off")
plt.show()
img_hyb = hybrid_image(img1, img2, 4, 2)
plt.figure(figsize=(15, 10))
plt.title("Highpass Plane and Lowpass Bird")
plt.imshow(vis_hybrid_image(img_hyb))
plt.axis("off")
plt.show()
img_hyb = hybrid_image(img2, img1, 5, 1)
plt.figure(figsize=(15, 10))
plt.title("Highpass Bird and Lowpass Plane")
plt.imshow(vis_hybrid_image(img_hyb))
plt.axis("off")
plt.show()
Cat and Dog
img1 = load_image("data/ex03/cat.bmp")
img2 = load_image("data/ex03/dog.bmp")
plt.figure(figsize=(15, 10))
plt.title("Cat and Dog")
plt.imshow(np.hstack([img1, img2]))
plt.axis("off")
plt.show()
img_hyb = hybrid_image(img1, img2, 3, 4)
plt.figure(figsize=(15, 10))
plt.title("Highpass Dog and Lowpass Cat")
plt.imshow(vis_hybrid_image(img_hyb))
plt.axis("off")
plt.show()
img_hyb = hybrid_image(img2, img1, 1.5, 2)
plt.figure(figsize=(15, 10))
plt.title("Highpass Cat and Lowpass Dog")
plt.imshow(vis_hybrid_image(img_hyb))
plt.axis("off")
plt.show()
Einstein and Marilyn
img1 = load_image("data/ex04/einstein.bmp")
img2 = load_image("data/ex04/marilyn.bmp")
plt.figure(figsize=(15, 10))
plt.title("Einstein and Marilyn")
plt.imshow(np.hstack([img1, img2]))
plt.axis("off")
plt.show()
img_hyb = hybrid_image(img1, img2, 1, 1)
plt.figure(figsize=(15, 10))
plt.title("Highpass Marilyn and Lowpass Einstein")
plt.imshow(vis_hybrid_image(img_hyb))
plt.axis("off")
plt.show()
img_hyb = hybrid_image(img2, img1, 2, 1.5)
plt.figure(figsize=(15, 10))
plt.title("Highpass Einstein and Lowpass Marilyn")
plt.imshow(vis_hybrid_image(img_hyb))
plt.axis("off")
plt.show()
Fish and Submarine
img1 = load_image("data/ex05/fish.bmp")
img2 = load_image("data/ex05/submarine.bmp")
plt.figure(figsize=(15, 10))
plt.figure(figsize=(15, 10))
plt.title("Fish and Submarine")
plt.imshow(np.hstack([img1, img2]))
plt.axis("off")
plt.show()
img_hyb = hybrid_image(img1, img2, 3, 2)
plt.figure(figsize=(15, 10))
plt.title("Highpass Submarine and Lowpass Fish")
plt.imshow(vis_hybrid_image(img_hyb))
plt.axis("off")
plt.show()
img_hyb = hybrid_image(img2, img1, 2, 1)
plt.figure(figsize=(15, 10))
plt.title("Highpass Fish and Lowpass Submarine")
plt.imshow(vis_hybrid_image(img_hyb))
plt.axis("off")
plt.show()
Durian and Ukulele
img1 = load_image("data/ex06/durian.jpg")
img2 = load_image("data/ex06/ukulele.jpg")
# Preprocessing the ukulele to create hybrid image
# Rotating the image and resizing to durian shape
img2 = cv2.resize(
nd.rotate(img2, -120, mode="constant", cval=1)[25:-60, :-25, :],
(img1.shape[1], img1.shape[0]),
)
# Clipping the outputs
img2[img2 < 0] = 0
img2[img2 > 1] = 1
plt.figure(figsize=(15, 10))
plt.imshow(img2.astype(np.uint8))
plt.title("Durian and Ukulele")
plt.imshow(np.hstack([img1, img2]))
plt.axis("off")
plt.show()
img_hyb = hybrid_image(img1, img2, 13, 3)
plt.figure(figsize=(15, 10))
plt.title("Highpass Ukulele and Lowpass Durian")
plt.imshow(vis_hybrid_image(img_hyb))
plt.axis("off")
plt.show()
img_hyb = hybrid_image(img2, img1, 20, 3)
plt.figure(figsize=(15, 7))
plt.title("Highpass Durian and Lowpass Ukulele")
plt.imshow(vis_hybrid_image(img_hyb))
plt.axis("off")
plt.show()
img1 = load_image("data/ex07/einstein.jpg")[10:-40, 25:-25, :]
img2 = load_image("data/ex07/marilyn.jpg")[:-35, :, :]
plt.figure(figsize=(15, 10))
# Reshaping the image to match each other
img1 = cv2.resize(img1, (img2.shape[1], img2.shape[0]))
plt.title("Einstein and Monroe Pt. 2")
plt.imshow(np.hstack([img1, img2]))
plt.show()
img_hyb = hybrid_image(img1, img2, 2, 1)
plt.figure(figsize=(15, 10))
plt.title("Highpass Marilyn and Lowpass Einstein")
plt.imshow(vis_hybrid_image(img_hyb))
plt.axis("off")
plt.show()
img_hyb = hybrid_image(img2, img1, 2, 1.5)
plt.figure(figsize=(15, 10))
plt.title("Highpass Einstein and Lowpass Marilyn")
plt.imshow(vis_hybrid_image(img_hyb))
plt.axis("off")
plt.show()